#include <iostream>
#include <fstream>
#include <vector>
#include <unistd.h>
#include <iostream>

#include "Vector.h"
#include "Color.h"

#define SQR(x) ((x)*(x))
#define _B "\033[1;38m"
#define B_ "\033[0m"

using namespace std;


class Particle {
	public:
		
		double m;
		Vector<double> R;
		double r;
		int id;
		
		Particle(): r(0) {}
		
		void print() {
			printf("P::print()\n");
			printf(" ├ m = %g\n", m);
			printf(" └ r = { %g, %g, %g }\n", R[0], R[1], R[2] );
		}
};


class Cell {
public:
	double m;
	Vector<double> R0;
	Vector<double> C;
	Vector<double> R1;
	
	vector<Particle*> p;
	int id;

	Cell():id(-2) { }
};


class PoissonNoise {
	public:
		Vector<double> R0;
		Vector<double> R1;
		
		PoissonNoise() :
		R0( Vector<double>(-60,-35,0) ),
		R1( Vector<double>(+60,+35,0) )
		{}
		
		template <typename T> 
		void set_box_r(vector<T> &p){
			
			Vector<> DR = R1-R0; 
			for (int i=0; i<p.size(); i++)
				p[i].R = R0 + Vector<>(	DR.x*drand48(), 
												DR.y*drand48(), 
												DR.z*drand48()	);
		}
};




class FlatGrid {
public:
	
	int nops;
	
	Vector<int> n;
	Vector<double> R0, R1;
	
	FILE *gp;
	
	FlatGrid():
	R0( Vector<double>(-60,-35,0) ),
	R1( Vector<double>(+60,+35,0) ),
 	n( Vector<int>(5*120,5*70,1) ),
 	nops(0),
	gp(popen("gnuplot --persist", "w"))
	{ }
	
	void create() {
		
		cell.resize(n.x*n.y*n.z);
		
		Vector<double> DR((R1.x-R0.x)/n.x, 
								(R1.y-R0.y)/n.y, 
								(R1.z-R0.z)/n.z);
		
		for (int k=0; k<n.z; k++)
		for (int j=0; j<n.y; j++)
		for (int i=0; i<n.x; i++) {
			
			Cell &cellijk = cell[i + n.x*j + n.x*n.y*k];
			cellijk.R0 = R0+Vector<double>(i*DR.x, j*DR.y, k*DR.z);
			cellijk.C  = R0+Vector<double>((i+0.5)*DR.x, (j+0.5)*DR.y, (k+0.5)*DR.z);
			cellijk.R1 = R0+Vector<double>((i+1)*DR.x, (j+1)*DR.y, (k+1)*DR.z);
			//cellijk.C = cellijk.R0 + 0.5*(cellijk.R1-cellijk.R0);
		}
	
	}
	
	template <typename T> 
	void map(vector<T> &p){
		
		Vector<double> DR(R1.x-R0.x, R1.y-R0.y, R1.z-R0.z);
		Vector<double> iDR( DR.x!=0 ? n.x/DR.x : 0, 
								  DR.y!=0 ? n.y/DR.y : 0,
								  DR.z!=0 ? n.z/DR.z : 0 );
		
		for (int i=0; i<p.size(); i++) {
			
			Vector<double> R = p[i].R - R0;
			
			int ix = R.x*iDR.x;
			int iy = R.y*iDR.y;
			int iz = R.z*iDR.z;
			
			int ijk = ix + n.x*iy + n.x*n.y*iz;
			Cell &cellijk = cell[ijk];
			
			p[i].id = ijk;
			cellijk.id = i;
			cellijk.p.push_back(&p[i]);
			
		}
	}
	
	vector<Cell*> cellact;
	template <typename T> 
	void init_rad_transfer(vector<T> &p) {
		
		map(p);
		
		// active cells
		for (int i=0; i<p.size(); i++)
			cellact.push_back(&cell[p[i].id]);
	}
	
	template <typename T> 
	void iterate_rad_transfer(vector<T> &p) {
		
		double cellsizex = (R1.x - R0.x)/n.x;
		double cellsizey = (R1.y - R0.y)/n.y;
		
		vector<Cell*> cellact_new;
		cellact_new.resize(0);
		
		for (int i=0; i<p.size(); i++)
			p[i].r+=cellsizex;
			
		
		
		for (int i=0; i<cellact.size(); i++) {
			
			bool freeCellsF = false;
			
			Cell &cellacti = *cellact[i];
			Particle &pi = p[cellacti.id];
		
			int ix = (cellacti.C.x-R0.x)/cellsizex;
			int iy = (cellacti.C.y-R0.y)/cellsizey;
			
			int ixx0 = ix-1; if (ixx0<0) ixx0 = 0;
			int ixx1 = ix+1; if (n.x<=ixx1) ixx1 = n.x-1;
			int iyy0 = iy-1; if (iyy0<0) iyy0 = 0;
			int iyy1 = iy+1; if (n.y<=iyy1) iyy1 = n.y-1;
			
			for (int iyy=iyy0; iyy<=iyy1; iyy++) 
			for (int ixx=ixx0; ixx<=ixx1; ixx++)
			if ( ixx!=ix || iyy!=iy) {
				
					
				int ixxyy = ixx + n.x*iyy;
				Cell &cellixxyy = cell[ixxyy];
				
				// 1. neiboring cell must be not ocupied 
				if (cellixxyy.id==-2) {
					// 2. neiboring cell must be inside expanding circle
					if ( (pi.R-cellixxyy.C).sqr()<=SQR(pi.r) ) {
						
						cellixxyy.id=cellacti.id;
						cellact_new.push_back(&cellixxyy);
						
					} else freeCellsF = true;
					nops++;
					
				} else if (cellixxyy.id!=cellacti.id) { // if we meet interface we make correction
				
						// who is closer to cellixxyy  p[cellacti.id] or p[cellixxyy.id] ?
						double dr2p0 = (cellixxyy.C - p[cellacti.id].R).sqr();
						double dr2p1 = (cellixxyy.C - p[cellixxyy.id].R).sqr();
					
						if (dr2p0<dr2p1)
							cellixxyy.id=cellacti.id;	
						else 
							cellixxyy.id=cellixxyy.id;	
						
						nops++;
					}
				
			}
			if (freeCellsF) cellact_new.push_back(&cellacti);
		}
		
		cellact.resize(0);
		cellact = cellact_new;
		
	}
	
	
	vector< Cell > cell;
	
	template <typename T> 	
	void plot(vector<T> &p, const char *fname) {
	
		double W = 1920;
		double H = 1080;
		
		
		fprintf(gp, "set xrange [%g:%g]\n", R0.x, R1.x);
		fprintf(gp, "set yrange [%g:%g]\n", R0.y, R1.y);

		
		
		double rm = 10;
		double lm = 10;
		double tm = 1;
		double bm = 1;
		
		double hr = 40;
		double wr = hr*(W-lm-rm)/(H-tm-bm);
		
		double boxsize = 0.10*(H-tm-bm)/n.y; // 1 - 0.165 2 - 0.15
		
		
		
		//fprintf(gp, "set term wxt size %g, %g\n", W, H);
		fprintf(gp, "set term pngcairo size %g, %g enhanced\n", W, H);
		fprintf(gp, "set output '%s'\n", fname);
		
		
		//fprintf(gp, "set lmargin at screen %g\n", lm/W) ;
		//fprintf(gp, "set rmargin at screen %g\n", (W-rm)/W);
		
		fprintf(gp, "set tmargin at screen %g\n", tm/H);
		fprintf(gp, "set bmargin at screen %g\n", (H-tm)/H);
		
		//fprintf(gp, "set border lc -1 lw 2\n");
		
		fprintf(gp, "set xtics %g, %g\n", R0.x, R1.x-R0.x);
		fprintf(gp, "set ytics %g, %g\n", R0.y, R1.y-R0.y);
		
		fprintf(gp, "set mxtics %i\n", n.x);
		fprintf(gp, "set mytics %i\n", n.y);
		
		fprintf(gp, "set format x ''\n");
		fprintf(gp, "set format y ''\n");
		fprintf(gp, "set object 1 rectangle from screen 0,0 to screen 1,1 fillcolor rgb '#FF5555' back\n");
		
		//fprintf(gp, "set grid mxtic mytics back lt 1, lt 3 lc rgb '#444444'\n");
		
		//fprintf(gp, "set size ratio -1\n");
		
		fprintf(gp, "set label 100000 \"nops =%7i\\nncells =%7i\\nnpoints=%7i\" at screen 0.05, 0.95 front font 'Arial, 18'\n", nops, cell.size(), p.size());
		
		
		/*
		for (int i=0; i<cell.size(); i++) //if (cell[i].id!=-2)
		//fprintf(gp, "set object %i rect from %g,%g to %g,%g fs solid 1.0 noborder fc lt %i back\n", i+1, cell[i].R0.x, cell[i].R0.y, cell[i].R1.x, cell[i].R1.y, cell[i].id );
		fprintf(gp, "set object %i rect from %g,%g to %g,%g fs solid 1.0 border lc rgb '#AAAAAA' lw 0.5 fc lt %i back\n", i+1, cell[i].R0.x, cell[i].R0.y, cell[i].R1.x, cell[i].R1.y, cell[i].id );
		*/
		
		/*
		for (int i=0; i<cell.size(); i++)
		fprintf(gp, "unset object %i \n", cell.size()+i+1 ); 
		
		
		for (int i=0; i<cellact.size(); i++)
		fprintf(gp, "set object %i rect from %g,%g to %g,%g fs solid 1.0 border lc rgb '#AAAAAA' lw 0.5 fc lt %i back\n", cell.size()+i+1, cellact[i]->R0.x, cellact[i]->R0.y, cellact[i]->R1.x, cellact[i]->R1.y, 1 );	
		*/
		/*
		for (int i=0; i<cell.size(); i++)
			for (int j=0; j<cell[i].p.size(); j++)
				fprintf(gp, "set object %i circle at %g,%g size %g fc rgb 'navy' \n", 2*cell.size()+i+1, cell[i].p[j]->R.x, cell[i].p[j]->R.y, cell[i].p[j]->r );
		*/
		
		fprintf(gp, "plot\ ");
		fprintf(gp, "'-' u 1:2:3 w p pt 5 ps %g lc rgb var t '',\ ", boxsize);
		fprintf(gp, "'-' u 1:2:3 w p pt 5 ps %g lc rgb var t '',\ ", boxsize);
		fprintf(gp, "'-' u 1:2:3 w p pt 7 ps 0.5 lc rgb var t ''\n");

		int vn = 2;
		int sn = (double)p.size()/vn;
		
		// all grid cells
		for (int i=0; i<cell.size(); i++) if (cell[i].id!=-2) {
			
			int vid = (double)cell[i].id/sn;
			int sid = cell[i].id - vid*sn;
			
			Color color; color.HSV2RGB(1.0*double(sid)/sn, 1, 1.0-0.4*double(vid)/vn);
			//Color color; color.HSV2RGB(0.0+1.0*double(cell[i].id)/(p.size()-1), 1.0, 1 );
			
			if (cell[i].id==-2) color.HSV2RGB(0, 0, 0.2);
			fprintf(gp, "%e %e %i\n", cell[i].C.x, cell[i].C.y, color.getColor() );
		}
		fprintf(gp, "e\n");
		
		// active cells 
		for (int i=0; i<cellact.size(); i++) {

			int vid = (double)cellact[i]->id/sn;
			int sid = cellact[i]->id - vid*sn;
			
			Color color; color.HSV2RGB(1.0*double(sid)/sn, 0.4-0.3*double(vid)/vn, 1 );
			//Color color; color.HSV2RGB(0.0+1.0*double(cellact[i]->id)/p.size(), 0.5, 1 );
			
			fprintf(gp, "%e %e %i\n", cellact[i]->C.x, cellact[i]->C.y, color.getColor());
		}
		fprintf(gp, "e\n");
		
		// mesh generation points
		for (int i=0; i<cell.size(); i++)
			for (int j=0; j<cell[i].p.size(); j++) {
				
				int vid = (double)cell[i].id/sn;
				int sid = cell[i].id - vid*sn;
			
				Color color; color.HSV2RGB(1.0*double(sid)/sn, 1.0, 0.5-0.5*double(vid)/vn);
				
				//Color color; color.HSV2RGB(0.0+1.0*double(cell[i].id)/p.size(), 0.6, 0.5 );
				fprintf(gp, "%e %e %i\n", cell[i].p[j]->R.x, cell[i].p[j]->R.y, color.getColor());
			}
		fprintf(gp, "e\n");
		
		
		
	}
	
	void print() {
		printf("\n");
		printf(_B"FlatGrid::print()\n"B_);
		printf(" ├ cellsize = {%20.16e, %20.16e, %g}\n", (R1.x - R0.x)/n.x, (R1.y - R0.y)/n.y, 0);
		printf(" ├ R0 = {% 20.16e, % 20.16e, %g}\n", R0.x, R0.y, R0.z);
		printf(" └ R1 = {% 20.16e, % 20.16e, %g}\n", R1.x, R1.y, R1.z);
		printf("\n");
	}
	
	
};



int main(int argc, char *argv[]) {
	
	//srand48(42);
	vector<Particle> p(90);
	
	PoissonNoise noise;
	noise.set_box_r(p);
	
	
	for (int i=0; i<p.size(); i++) p[i].R.z=0;
	
	FlatGrid flatgrid;
	
	flatgrid.create();
	flatgrid.init_rad_transfer(p);
	flatgrid.print();
	
	char fname[256];
	
	int k = 0, kk = 0;
	for (int i=0; i<300; i++) {
		
		for (int j=0; j< (i==0 ? 25 : 1); j++) {	
			printf("%04i / %04i\n", i, k);
			sprintf(fname, "/home/yurinds/tmp/voronoi/%04d.png", k);
			flatgrid.plot(p, fname);		
			k++;
			if (flatgrid.cellact.size()==0) kk++;
		}
		
		flatgrid.iterate_rad_transfer(p);
		
		
		if (100<kk) break;
		
	}
	
	return 0;
	
}